home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_01 / bnarth3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-17  |  5.4 KB  |  258 lines
  1. /*
  2.  *   MIRACL arithmetic routines 3 - powers and roots
  3.  *   bnarth3.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include "miracl.h"
  8. #define abs(x)  ((x)<0? (-(x)) : (x))
  9.  
  10. /* Access global variables */
  11.  
  12. extern int   depth;    /* error tracing ..*/
  13. extern int   trace[];  /* .. mechanism    */
  14. extern small base;     /* number base     */
  15. extern int   lg2b;     /* bits in base    */
  16. extern small base2;    /* 2^lg2b          */
  17. extern bool  check;    /* overflow check  */
  18.  
  19. extern big w0;
  20. extern big w1;        /* workspace variables */
  21. extern big w2;
  22. extern big w3;
  23. extern big w4;
  24. extern big w5;
  25. extern big w6;
  26.  
  27. int logb2(x)
  28. big x;
  29. { /* returns number of bits in x */
  30.     int xl,lg2;
  31.     if (ERNUM || size(x)==0) return 0;
  32.     depth++;
  33.     trace[depth]=49;
  34.     if (TRACER) track();
  35.     copy(x,w0);
  36.     insign(PLUS,w0);
  37.     lg2=0;
  38.     xl=w0[0];
  39.     if (base==base2)
  40.     {
  41.         lg2=lg2b*(xl-1);
  42.         shift(w0,1-xl,w0);
  43.     }
  44.     else while (w0[0]>1)
  45.     {
  46.         subdiv(w0,base2,w0);
  47.         lg2+=lg2b;
  48.     }
  49.     while (w0[1]>1)
  50.     {
  51.         lg2++;
  52.         w0[1]/=2;
  53.     }
  54.     depth--;
  55.     return lg2;
  56. }
  57.  
  58. void expb2(x,n)
  59. big x;
  60. int n;
  61. { /* sets x=2^n */
  62.     int i,r;
  63.     if (ERNUM) return;
  64.     convert(1,x);
  65.     if (n==0) return;
  66.     depth++;
  67.     trace[depth]=50;
  68.     if (TRACER) track();
  69.     if (n<0)
  70.     {
  71.         berror(10);
  72.         depth--;
  73.         return;
  74.     }
  75.     r=n/lg2b;
  76.     if (base==base2)
  77.     {
  78.         shift(x,r,x);
  79.         x[x[0]]<<=(n%lg2b);
  80.     }
  81.     else
  82.     {
  83.         for (i=1;i<=r;i++)
  84.             premult(x,base2,x);
  85.         premult(x,((small)1<<(n%lg2b)),x);
  86.     }
  87.     depth--;
  88. }   
  89.  
  90. void power(x,n,z,w)
  91. big w;
  92. big x;
  93. int n;
  94. big z;
  95. { /* raise big number to int power  w=x^n *
  96.    * (mod z if z and w distinct)          */
  97.     copy(x,w5);
  98.     zero(w);
  99.     if(ERNUM || size(w5)==0) return;
  100.     convert(1,w);
  101.     if (n==0) return;
  102.     depth++;
  103.     trace[depth]=17;
  104.     if (TRACER) track();
  105.     if (n<0)
  106.     {
  107.         berror(10);
  108.         depth--;
  109.         return;
  110.     }
  111.     if (w!=z) divide(w5,z,z);
  112.     if (w==z) forever
  113.     { /* "Russian peasant" exponentiation */
  114.         if (n%2!=0) multiply(w,w5,w);
  115.         n/=2;
  116.         if (ERNUM || n==0) break;
  117.         multiply(w5,w5,w5);
  118.     }
  119.     else forever
  120.     {
  121.         if (n%2!=0) mad(w,w5,w5,z,z,w);
  122.         n/=2;
  123.         if (ERNUM || n==0) break;
  124.         mad(w5,w5,w5,z,z,w5);
  125.     }
  126.     depth--;
  127.     return;
  128. }
  129.  
  130. void powmod(x,y,z,w)
  131. big w;
  132. big x;
  133. big y;
  134. big z;
  135. {  /*  calculates w=x^y mod z */
  136.     if (ERNUM) return;
  137.     copy(y,w1);
  138.     copy(x,w2);
  139.     zero(w);
  140.     if (size(w2)==0) return;
  141.     convert (1,w);
  142.     if (size(w1)==0) return;
  143.     depth++;
  144.     trace[depth]=18;
  145.     if (TRACER) track();
  146.     if (size(w1)<0) berror(10);
  147.     if (w==z)       berror(7) ;
  148.     if (ERNUM)
  149.     {
  150.         depth--;
  151.         return;
  152.     }
  153.     divide(w2,z,z);
  154.     forever
  155.     { /* "Russian peasant" exponentiation */
  156.         if(subdiv(w1,2,w1)!=0)
  157.             mad(w,w2,w2,z,z,w);
  158.         if(ERNUM || size(w1)==0) break;
  159.         mad(w2,w2,w2,z,z,w2);
  160.     }
  161.     depth--;
  162. }
  163.  
  164. bool root(x,n,w)
  165. big x,w;
  166. int n;
  167. {  /*  extract  lower approximation to nth root   *
  168.     *  w=x^(1/n) returns TRUE for exact root      *
  169.     *  using Newtons method                       */
  170.     int sx,dif,s,p,d,lg2,lgx,rem;
  171.     bool full;
  172.     if (ERNUM) return FALSE;
  173.     if (size(x)==0 || n==1)
  174.     {
  175.         copy(x,w);
  176.         return TRUE;
  177.     }
  178.     depth++;
  179.     trace[depth]=16;
  180.     if (TRACER) track();
  181.     if (n<1) berror(11);
  182.     sx=exsign(x);
  183.     if (n%2==0 && sx==MINUS) berror(9);
  184.     if (ERNUM) 
  185.     {
  186.         depth--;
  187.         return FALSE;
  188.     }
  189.     insign(PLUS,x);
  190.     lgx=logb2(x);
  191.     if (n>lgx)
  192.     { /* root must be 1 */
  193.         insign(sx,x);
  194.         convert(sx,w);
  195.         depth--;
  196.         if (lgx==0) return TRUE;
  197.         else        return FALSE;
  198.     }
  199.     expb2(w2,1+lgx/n);           /* guess root as 2^(log2(x)/n) */
  200.     s=(-((x[0]-1)/n)*n);
  201.     shift(w2,s/n,w2);
  202.     lg2=logb2(w2);
  203.     full=FALSE;
  204.     if (s==0) full=TRUE;
  205.     d=0;
  206.     p=1;
  207.     while (!ERNUM)
  208.     { /* Newtons method */
  209.         copy(w2,w3);
  210.         shift(x,s,w4);
  211.         check=OFF;
  212.         power(w2,n-1,w6,w6);
  213.         check=ON;
  214.         divide(w4,w6,w2);
  215.         rem=size(w4);
  216.         subtract(w2,w3,w2);
  217.         dif=size(w2);
  218.         subdiv(w2,n,w2);
  219.         add(w2,w3,w2);
  220.         p*=2;
  221.         if(p<lg2+d*lg2b) continue;
  222.         if (full && abs(dif)<n)
  223.         { /* test for finished */
  224.             while (dif<0)
  225.             {
  226.                 rem=0;
  227.                 decr(w2,1,w2);
  228.                 check=OFF;
  229.                 power(w2,n,w6,w6);
  230.                 check=ON;
  231.                 dif=compare(x,w6);
  232.             }
  233.             copy(w2,w);
  234.             insign(sx,w);
  235.             insign(sx,x);
  236.             depth--;
  237.             if (rem==0 && dif==0) return TRUE;
  238.             else                  return FALSE;
  239.         }
  240.         else
  241.         { /* adjust precision */
  242.             d*=2;
  243.             if (d==0) d=1;
  244.             s+=d*n;
  245.             if (s>=0)
  246.             {
  247.                 d-=s/n;
  248.                 s=0;
  249.                 full=TRUE;
  250.             }
  251.             shift(w2,d,w2);
  252.         }
  253.         p/=2;
  254.     }
  255.     depth--;
  256.     return FALSE;
  257. }
  258.